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Abstract 

The Laplacian growth (the Hele-Shaw problem) of multiply-connected domains in 
the case of zero surface tension is proven to be equivalent to an integrable system 
of Whitham equations known in soliton theory. The Whitham equations describe 
slowly modulated periodic solutions of integrable hierarchies of nonlinear differential 
equations. Through this connection the Laplacian growth is understood as a flow in 
the moduli space of Riemann surfaces. 



*Department of Mathematics, Columbia University, New York, USA, Landau Institute and ITEP, 
Moscow, Russia 

tLos Alamos National Laboratory, MS-P365, Los Alamos, NM 87545, USA 

■fjames Frank Institute and Enrico Fermi Institute of the University of Chicago, 5640 S.Ellis Avenue, 
Chicago, IL 60637, USA and Landau Institute for Theoretical Physics, Moscow, Russia 

§CNLS, LANL, Los Alamos, NM 87545, USA, and Institute of Biochemical Physics, 4 Kosygina st., 
119991, Moscow, Russia and ITEP, 25 B.Chercmushkinskaya, 117259, Moscow, Russia 



1 



1 Introduction 



It is not uncommon that nonlinear differential equations which possess an intcgrablc struc- 
ture emerge in important problems of hydrodynamics [1]. The Korteweg de Vrics equation 
describing nonlinear waves in dispersive media is perhaps the most familiar example. 

In recent years, integrable structures were found in another class of hydrodynamics prob- 
lems leading to a pattern formation in a regime far from equilibrium [2]. Growth problems 
of this type are unified by the name Laplacian growth. In this paper, we further develop a 
hnk between growth processes and sohton theory. We extend the results of Ref.[2] to the 
case of multiply-connected domains and identify the set of growth processes with a universal 
Whitham hierarchy of integrable equations. The latter unveils the mathematical structure 
of the growth and set a place for growth models in the realm of soliton theory. 

Laplacian growth, also known as the Hele-Shaw problem, refers to dynamics of a moving 
front (an interface) between two distinct phases driven by a harmonic scalar field. This field 
is a potential for the growth velocity field. The Laplacian growth problem appears in different 
physical and mathematical contexts and has a number of important practical applications. 
The most known ones are filtration processes in porous media, viscous fingering in the Hele- 
Shaw cell, electrodeposition and solidification in undercooled liquids. A comprehensive list 
of relevant papers published prior to 1998 can be found in [3]. 

The most interesting and the most studied dynamics occurs in two-dimensional spatial 
geometry. To be definite, we shall speak about an interface between two incompressible 
fiuids with very different viscosities on the plane. In practice the 2D geometry is realized in 
a Hele-Shaw cell - a narrow gap between two parallel plates (Fig. 1). In this version, the 
problem is also known as the Saffman- Taylor problem or viscous fingering. For a review, see 
[4]. Importance of studies of the Laplacian growth with more than one nonviscous droplet 
speaks for itself. When rates of flow are considerable, then, because of fingering instability, 
new droplets are pinched off and change their shapes, so the whole dynamics considerably 
changes in comparison with a single bubble dynamics. (See, e.g., experimental works [5].) 

To be more precise, consider the case when there are several disconnected domains in the 
Hele-Shaw cell occupied by a fluid with low viscosity (water). We call them water droplets. 
Their exterior, which is in general a multiply-connected domain, is occupied by a viscous 
fluid (oil). All components of the oil/water interface are assumed to be smooth curves. Oil 
is sucked out with flxed rates Qj through sinks placed at some flnite points aj or at inflnity 
(edges of the Hele-Shaw cell). Water is injected into each water droplet with rates qa, some 
of which may be negative or equal to zero. Hele-Shaw flows in general setting, including the 
multiply-connected case, were discussed in [6]-[ll]. 

Viscous flows are governed by gradient of the pressure field in the fluids. In the oil 
domain, the local velocity V = {Vx, Vy) of the fluid is proportional to the gradient of pressure 

— * 

p = p{x,y) (Darcy's law): V — —nVp, where k is called the filtration coefficient. In what 
follows, we choose units in such a way that n — ^. In particular, the Darcy law holds on the 
outer side of the interface thus governing its dynamics: 
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Figure 1: The Hele-Shaw cell (arrows show direction of the Bow forced by a pump). 



Here dn is the normal derivative. 

This simple dynamics results to complicated unstable patterns often growing beyond 
control. The most recent experimentally produced pattern can be seen in Ref. [5]. 

In this paper we discuss only the idealized problem, namely that surface tension equals 
zero. We call it the idealized Laplacian growth (ILG). Zero surface tension means that 
pressure does not jump across the boundary. Assuming that viscosity of water is small 
enough comparing to the viscosity of oil, pressure is constant inside each water droplet. 
However, the values of pressure may be different in different droplets and may also depend 
on time. Let be pressure in the a-th droplet, then zero surface tension implies that p — Pa 
on the outer side of the interface as well. 

Since fluids are incompressible {VV = 0) the Darcy law implies that the pressure field p 
is a harmonic function in the exterior (oil) domain except at the points where the oil pumps 
are located. In the case of zero surface tension, pressure is a solution of the time-dependent 
boundary problem for the Laplace equation with p = pa on the boundary components. The 
interface moves according to the Darcy law (1), so that the boundary problem changes with 
time. Note that the problem is non-local since the gradient of pressure around boundary 
depends on the shape of the domain whole. 

When the interface bounds a simply-connected domain, an effective tool for dealing with 
the Laplacian growth is the time dependent conformal mapping technique (see e.g. [4]). 
Passing to the complex coordinates z — x + iy, z — x — iy on the physical plane, one may 
describe dynamics in terms of a moving conformal map from a simple reference domain, say 
the unit disk in a "mathematical plane" , onto a growing domain in the physical plane. 

If the interface has several disconnected components the conformal map approach meets 
fundamental difficulties. Uniformizing maps of multiply-connected domains are essentially 
more complicated mainly because there are no simple reference domains and, moreover, any 
possible reference domain depends on the growing domain itself. 

Instead, given growing domains, one may construct an algebraic curve over complex 
numbers, or a Riemann surface, endowed with some additional structures. Then the Darcy 
law translates the growth to an evolution of the algebraic curve, thus providing a remarkable 
family of flows on the (extended) moduli space of Riemann surfaces. 

This is the meeting point of the growth processes and soliton theory. Wc show that the 
evolution of the Riemann surface built upon the moving interface in the Hele-Shaw cell is 
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governed by the universal Whitham hierarchy of sohton theory. It generahzes the disper- 
sionless 2D Toda hierarchy which was shown in Ref. [2] to describe the moving boundary 

problems in the simply-connected case. The Whitham equations have been originally intro- 
duced to describe slow modulations of periodic solutions to differential equations [12]. Their 
relation to the algebraic geometry, using example of the KdV equation, was found in [13]. 
In general setting, the universal Whitham hierarchy was introduced in [14, 15, 16]. 

The mathematical equivalence between the Laplacian growth with zero surface tension 
and the Whitham equations, established in this work, allows one to treat the former as a 
singular limit of dispersive waves obeying soliton equations. Such a singular limit may often 
lead to solutions which develop singularities within finite time. 

In the context of the Laplacian growth, the singularities show up as cusps generated 
by an initially smooth interface [19], after which the idealized description no longer makes 
any physical sense. This feature signifies that the very problem with zero surface tension is 
ill-posed. 

A similar difficulty has been known to exist for the Whitham equations (see, e.g., [17]). 
Some of their solutions, being initially regular, sooner or later become singular and so can not 
be extended to all times. The Korteweg de Vries equation with zero dispersion (the Hopf 
equation) is the most familiar example. Here any smooth decreasing Cauchy data evolve 
into a "shock wave" with an overturned front which is physically meaningless (a "gradient 
catastrophe" ) . This simply means that the dispersionless approximation does not work in a 
vicinity of the catastrophe. Similarly, the zero surface tension assumption is not valid in a 
vicinity of the cusp formation. 

These singularities are in fact artificial and can be successfully resolved by methods 
developed in the theory of slow modulations of exact periodic solutions to soliton equations 
[1, 18]. In subsequent works we hope to apply these methods to the Laplacian growth using 
the proven below equivalence between the two disciplines. 

2 Linearization of the ILG dynamics 

Remarkably, the ILG dynamics, initially formulated as a non-local and highly nonlinear 
problem, admits an exact linearization in the space of harmonic moments of the viscous 
domain. By linearization we mean here a change of variables which converts the nontrivial 
LG dynamics into a simple linear one. A familiar (but rather loose) analogy is passing to 
action-angle variables in classical mechanics or the inverse scattering transform in the soliton 
theory. 

Our starting point is the fact that the ILG is a simple linear shift in the space of harmonic 
moments of the growing domain. This statement goes back to seminal Richardson's paper 
[6]. In that paper, it was shown that if there is the only sink at infinity, all moments are 
conserved except the moment of constant function (the area of the droplet) which changes 
linearly with time. In fact, it is absolutely clear that for incompressible fluids and fixed 
pumping rates areas of the droplets, if change at all, always do this linearly with time. For 
arbitrary location of the sink, and also for several sinks at different points, a simple extension 
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of this result states that in general all moments change linearly with time (with different 
coefficients which may be zero). 

In the case of several water droplets, the set of harmonic moments should be supplemented 
by a finite number of extra parameters, one for each extra droplet, which are basically mo- 
ments of harmonic functions with multivalued analytic parts. This set of variables is enough 
to characterize the geometry of the growing domain. Alternatively, the new parameters may 
be areas of the water droplets. Depending on which type of external physical conditions in 
the water droplets is realized (fixed pressure differences or fixed pumping rates), the ILG 
dynamics becomes linear either in the former or in the latter variables. 

Most of the material of Sees. 2.1 and 2.2 is spread through the literature (see, e.g., 
[6, 7, 8, 9] and Chapter 5 of the book [11]). To make our exposition self-contained, we review 
them from a unifying point of view. 

2.1 The time-dependent boundary problem 

Consider an ILG process with the point-like oil pumps with powers Qj at some points 
located far enough from the moving interface. Mathematically this means 



where Cj is a small contour encircling the point a^, ds stands for the differential of the arc 
length, Vri is the component of the fluid velocity normal to the contour, with the normal 
vector pointing outside the circle Cj. Both Qj and aj are assumed to be time- independent. 
Oil may be also sucked at infinity. Physically this means, for example, that oil is removed 
from the edge of a large Hele-Shaw cell. Mathematically one puts one of aj equal to infinity 
and defines the pumping rate at infinity, Qoo, as 



Here Cqo is a big contour encircling the whole system of water droplets and all the point-like 
pumps, if any. The oil pumping rates are assumed to be positive when oil is sucked and 
negative if it is injected into the Hclc-Shaw cell. 

One can also consider extended sources or sinks of oil, for instance, continuously dis- 
tributed along lines, like in Ref. [9]. To avoid irrelevant technical complications, we consider 
point-like oil pumps only, giving brief remarks on the more general case when necessary. 

Let Doii be the region of the plane occupied by oil (an infinite domain containing infinity), 
and D be the region occupied by water. We assume that there are g -\- 1 water droplets in 
the Hclc-Shaw cell, which arc compact domains bounded by smooth non- intersecting curves. 
Let Da be the a-th water droplet, a = 0,1, ... ,g, so that D is their union (Fig. 2). It is 
convenient to think of the 0-th droplet as of the "main" one, having in mind that it is this 
droplet that remains in the system in the simply-connected case. Let us assume, just for a 
notational convenience, that the origin is in the zero droplet, e Dq. 




(2) 
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Figure 2: Water droplets (g — 3). 



The incompressibility implies that the oil pumps are able to work only if there are some 
sources of water which supply water to at least one of the droplets. We call them water 
pumps and introduce the pumping rate Qa of water in the a-th droplet: 

f Vnds = TTQa , (3) 

Here the normal vector points outside D„. The pumping rate is assumed to be positive 
when water is injected into the a-droplet. Obviously, the oil and water pumping rates are 
constrained by the relation 

j " 

The pressure field in oil obeys the equation 

Ap{z)^A7rj:Q,6^'\z-a,), (4) 
i 

where 6^'^\z) is the two-dimensional delta-function and A = d"^ + dy = Ad^dz is the Laplace 
operator on the plane. In other words, p = p{z) is a harmonic function in Dqu with the 
asymptote 

p{z) — Qj log \z — Qjl"^ + . . . in the vicinity of Uj 

Taking into account the Darcy law, this agrees with the definition of the Qj (2). On the a-th 
component of the interface the pressure field takes a constant value Pa (which may depend 
on time). Without loss of generality, we set po = 0. 

To determine velocity of the interface, one should solve a time-dependent Dirichlet-like 
boundary problem. It has a unique solution, and so pressure in Doii is uniquely determined 
as soon as one specifies pressures Pa in water droplets. This means that with given p^s there 
is no freedom in water pumping. Indeed, the pumping rates in all water droplets are to be 
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determined from (3), which states that — — ^ §qq^ dnpds. li PaS are maintained constant 
with time, q^s are in general certain comphcated functions of time. Alternatively, one may 
control the pumping rates q^ keeping them constant, then pressures in water droplets are 
uniquely determined by the dynamical equations and in general exhibit a nontrivial time 

dependence. 

So the physical problem is not yet defined by the local growth law alone. One should 
add some physical conditions in water droplets. We distinguish two cases: 

(I) Fixed pressure differences Pa' Pa — 0, then q^ — qa{t) are to be determined; 
(II) Fixed pumping rates qa- qa = 0, then pa = Pa{t) are to be determined 

Here the dot means the time derivative. Various "mixed" conditions (say, when we fix 
some of PaS and some of q^s or linear combinations thereof) are not meaningless from the 
mathematical point of view but look somewhat artificial in the physical context. 

For an abuse of quantum-mechanical terminology, one may refer to the cases (I) and (II) 
as the LG problem in p-representation and ^-representation, respectively. 



2.2 Integral transforms of the oil domain 

The Cauchy transform. Here we closely follow Richardson's paper [8] . Let us consider 
the Cauchy transform of the oil domain: 

AW = V ^ (5) 

(A cut-off at some radius R is implied, at large R the integral does not depend on the cut-off 
which may be safely moved to infinity.) This function is continuous across all boundaries 
and analytic for z in each hole Dq, while in Doo the function h{z) — z is analytic. So we can 
write 

h'^{z) for z in Dq, 
hiz) = (6) 
^ z + h (z) for z in Doo 

where h^{z) is an analytic function in and h~{z) is analytic for z in Don- One may 
analytically continue the functions h" outside the regions where they are defined by the 
Cauchy transform. In general, all -|- 1 functions h'^iz) obtained in this way are expected 
to be different. 



The contour integral representation of these functions reads: 



2m JdD C — z \ r ■ r\ 

^ \ h [z) tor z m Don 

Note that h'^{z) — h~{z) = z on the a-th boundary curve. 



(7) 
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The time derivative h{z) = dth{z;t), is found straightforwardly using the integral repre- 
sentation (5) and the Darcy law: 

If 1 If 1 

h{z) = f Vnds ^ — f dnP{C)ds = 

where we subtracted and added the integral of pdn^^ over the whole boundary, and used 
the fact that p is constant along any component of the boundary. It is easy to sec that 
irrespectively of whether the point z is inside or outside D^, all contour integrals in the 
second term vanish. The first term can be transformed using the Green theorem: 

dnP - pdn^] ds^- j Ap(C) - p(C) A d\ 

Z-C Ju^At) \ z-C Z-C 



The Laplacian of is given by (4). If z is inside a water droplet, the function (regarded 
as a function of Q is harmonic for all C, in Doii and the second term vanishes. If z is in Doii, 



L 



Finally, we get: 



In the integrated form. 



p{Q^A—-\dX^^7:d,p{z) 
hl{z)^Y.-^ for^inD, (8) 

■ Qj Z 

h- {z) = V + d,p{z) for z in Don (9) 

ya,-z 

hUz;t) = hUz;0) + tY:-^ (10) 

■ Uj z 

We see that the increment of the is a rational function, and, moreover, it is the same 
for all a — 0, 1, . . . , g^. We also see that this function is entirely determined by the output 
powers and locations of oil pumps, no matter what conditions we impose on the pressures 
and pumping rates in water droplets. 

For a linearly distributed oil source, p{z) is a potential of a simple layer. Then one 
finds that the functions /i^, though no longer rational, are still analytic continuations of a 
single analytic function unless the support of the simple layer forms a noncontractable cycle 
encircling at least one of the water droplets. 

A remarkable property of the Cauchy transform is its linear dependence on time t (pro- 
vided the parameters of the oil pumps are time-independent). This is the key point that 
allows one to linearize the ILG dynamics. 
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The Coulomb potential. Along with the Cauchy transform of the domain Doii it is useful 
to consider the potential generated by fictitious 2D Coulomb charges uniformly distributed 
in Doii: 



/ log 











(11) 



(The same cut-off as in the Cauchy transform is implied.) Clearly, this function is harmonic 
in each water droplet and its 2;-derivative coincides with the h{z): dz(f){z) = h{z). Repeating 
the above calculation for (f){z) we get (for z in water droplets): 



1 







/ log 









1 - 



ds 



This yields: 



0(^) = Pa - Z] Qj log 



for z in 



(12) 



The Laplacian growth equation. For completeness, let us demonstrate that the Lapla- 
cian growth equation, usually derived, in the simply-connected case, using the time depen- 
dent conformal map technique, follows from the time derivative of the Cauchy transform. 

To this end, wc calculate the difference of the boundary values in two ways. On the 
one hand, it is obtained by subtracting (8) and (9): 

h'^{z) — h~{z) — —dzp{z) , for z on any boundary contour (13) 

On the other hand, 

' ^ ' 27li ' \JdD C-Zj 

can be found directly using a paramctrization of the family of contours z = z{(T,t), where a 
is a parameter along the contour. For each component of the boundary we have: 

Zdz \ _ r ( ~tZa + ZZ„t _ ZZ^Zt \ _ 

z — aj J \ z — a {z — aY j 

/. .1 \ f f ZtZcr — ZcrZf\ dz 

(mtegratmg by parts) = f — 

J \ z — a J Za 

The jump of the boundary values of the analytic function defined by the latter Cauchy 
integral is equal to — 'z^zt{zcr)~^ ■ Combining this with (13), we obtain the relation for 
differentials along the boundary curves, 

dtz{a, t)dz — dtz{a, t)dz — —dzp{z)dz (14) 

valid for any paramctrization of the contours. In a simply-connected case with an oil sink 
at infinity p{z) — — 2Q^ log \w{z)\, where w{z) is a conformal map of the oil domain to 
the exterior of the unit disk. Choosing a — —ilogw{z) Eq.(14) becomes the celebrated 
Laplacian growth equation 

It was first derived in [20]. 
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Figure 3: The domain Dqw is the exterior of the regions bounded by the dashed hnes. 

2.3 Dual systems of local coordinates in the space of multiply- 
connected domains 

Our next goal is to introduce special local coordinates in the space of multiply-connected 
domains, which evolve linearly in time. The time evolution of the Cauchy transform suggests 
that such coordinates are basically harmonic moments of the oil domain. 

The proper basis of harmonic functions. Let us consider a time-independent domain 
□oil C Doii with the same connectivity as Doq, as shown in Fig. 3, and define a proper basis 
of harmonic functions in Doq. A basis is said to be proper if any harmonic function in Doii 
is representable as a linear combination (possibly infinite) of the basis functions such that 
it converges everywhere in Doii- In the case of a single water droplet, the basis consisting of 
functions z~"' and their conjugates is clearly a proper one. However, this basis is no longer 
proper on the plane with more than one hole. Indeed, in this case one has to incorporate 
functions with singularities in any hole, not only in Dq, otherwise the series converges only 
in some simply-connected neighborhood of infinity. 

To construct a proper basis, we need some auxiliary data. Inside each water droplet, 
let us fix a point, Za G Dq, which does not move with time (Fig. 3). These points may be 
thought of as locations of point-like water pumps. Without any loss of generality, we set 
^0 = 0. 

Consider the following single- valued analytic functions in Doq: 




A; > 0, 



(15) 
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where the running index 7 is understood modulo g + 1, i.e., it is imphed that = 
At /c = we set -00 = 1- In a more exphcit form, we write 

0-1 

tp{g+l)m+a = Y[{z - ZpY^ , (16) 

/3=0 

where m > 0, a = 0, 1, . . . , gi (when a = the product is set to be 1) and 

r{z) = f{{z-zp) (17) 

13=0 

These functions constitute a proper basis of single- valued analytic functions. (As is seen from 
(16), the domain Dqu can be choosen to be \r{z)\ > C for some positive C.) This basis is a 
simplified version of the special Laurent-Fourier type basis on Riemann surfaces constructed 
in [21]. In the case of the single hole Dq with 2:0 = it coincides with the standard one: 

Any harmonic function is the real part of an analytic function. However, in the multiply- 
connected case these analytic functions are not necessarily single- valued (only their real parts 
have to be single- valued). The full basis of harmonic functions is, therefore: i^k{z), i^k{z), 
k > 'ipo — and g functions of the form 

2 



£a{z) = log 



Z 



« = (18) 



which have multivalued analytic part. 

An important feature of the basis ■0fe is the existence of the "dual" basis of differentials, 
dxk{z), such that 

TT-. f i^kdXn = (19) 

Explicitly, they are given by 

dz ^-'^ 

dxk{z) = 7 ^—TT ^dzl[{z- z^) (20) 

{z - Zk-l)ll^k{z) ^^io 

where is again understood as where k = {g + l)m + a with non-negative integer m and 
< a < (?. It is clear that x'ki^) — dxk/dz, for k >1, are polynomials of degree k — 1. 

One can easily check that the following expansion of the Cauchy kernel holds true: 

-1 00 

-r— = j:MOx'k{z) (21) 

^ k=i 

and the series converges if \r{z)\ < \r{C)\- In particular, the series certainly converges if ( is 
in Doii while z is close enough to any one of the points Za- A similar expansion holds for the 
logarithmic kernel: 

log(C - z) ^ log(C - Za) -J2 ^fc(C) / dxk for C in Dou and z in (22) 

k=i 

with the same condition of convergence. 
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The hcirmonic moments. Using the expansion of the Cauchy kernel (21), let us expand 
the function h'^{z) in a series, which converges at least for z close enough to any one of the 
points Za- We obtain: 

h^{z)^^knx',{z) 



k>l 



where 



(23) 



are normalized harmonic moments of the domain Doii- (In the case of the single water droplet 
Do they have the familiar form = —-^ . Their time derivatives are read from 

(8): " _^ " 

n^lEQjMaj) (24) 

In case of the only sink at infinity the r.h.s. vanishes [6]. The moment Tq is the area of Doii. 
It is infinite, but as long as we need only its change, we may equivalently consider the area 
of the complimentary domain 

To = - ^ Area (D,) 



TT 



a=0 



Evidently, Tq is equal to the total power of oil pumps: Tq = Qoo + J2j Qj- 

To characterize the geometry of a multiply-connected domain with g + 1 boundary com- 
ponents one needs g extra parameters. There are different possibilities to choose them. One 
of possible choices is as follows. Using the expansion of the logarithmic kernel (22), or inte- 
grating the expansion of the h{z), we can represent the Coulomb potential function in the 
form 



(l>{z) ^^^ + 2Re I Y^kTk T J 



for z in 



(25) 



where — 4>{za) arc integration constants. Since 0(0) = 0, we set $o = 0. The other inte- 
gration constants, . . . , may serve as the additional independent parameters. Clearly, 



= - / log 



2 



are harmonic moments of the Doii with respect to the functions la{z) (18). These moments 
are always real. The time derivative = 0(^a) is given by (12): 



Pa - X)Qi4(ai) 



(26) 



There is an alternative choice of the g extra parameters which is "dual" to the choice 
above. Let 

= (27) 

TT 
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be areas of the water droplets (divided by n) , then 



Sa — 



(28) 



by the definition of Qa- Clearly, J2a=o — Tq, so Si, . . . ,Sg can be taken as independent 
parameters. 

Local coordinates in the space of multiply-connected domains. The basic fact 
from the theory of deformations of planar domains (which we adopt without proof in this 
paper) is that the parameters T^, $q, or T^, can serve as local coordinates in the space of 
planar multiply-connected domains. This means that any deformation which preserves these 
parameters is trivial and any vector field in the space of these parameters generates a well- 
defined deformation of the initial domain (see [22] for details) . Recovery of the domain having 
these parameters constitutes the classical inverse potential problem [23] for the multiply- 
connected case. The fact that the coordinates introduced above are good ones amounts to 
local existence and uniqueness of such a domain. An effective reconstruction, however, is not 
generally feasible and even specific examples are challenging and merit separate attention. 

The formulas for time derivatives of T^, $q, and Sa tell us that any ILG fiow can be 
represented, in these coordinates, as a vector field with constant coefficients. 

To summarize, we have introduced two systems of local coordinates in the space of 
multiply-connected domains: 

• T^-coordinates: the harmonic moments Tq, Ti, ■ ■ ■ (and their complex conjugates) 
and $i,$2,...,^fl; 

• TjS-coordinates: the harmonic moments Tq, Ti, T2, . . . (and their complex conjugates) 
and Si, S2,..., Sg. 

Now, from (24), (26) and (28) it is clear that Richardson's result can be reformulated by 
saying that the ILG dynamics with fixed pressure differences (the p-representation) is lin- 
earized in the T^-coordinates while the ILG dynamics with fixed pumping rates (the q- 
representation) becomes linear in TjS'-coordinates. 

2.4 Elementary growth processes 

A linear time dependence of the local coordinates suggests to treat any general ILG process 
as a superposition of certain "elementary" processes. It is natural to associate with each 
elementary fiow its own time variable. 

In the p-representation, the elementary processes are: 

(pi) Oil is sucked from a point a with the unit rate {Q — 1), with maintaining equal 
pressures in all droplets {pa — 0): 



To — 1, Tk — T 



a 



4(a) 
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With this process we associate the time variable T^"^ such that 

(the values of the coefficients follow from (24), (26)). Note that T^"^ is the amount of 
oil sucked out from the point a during the process. 

(p2) Water is redistributed between the droplets by applying the unit pressure difference 
between the a-th and the 0-th droplets: pp — S^p, with no pumps in oil: 

To = Tk = 0, = 6a(3 

With this process we associate the time variable T^") such that 

d _ d 

It is the amount of water injected into the 0-th water droplet during the process. 
In the g-representation, the elementary processes are: 

(ql) Oil is sucked from a point a with the unit rate {Q = 1), with water being added to the 
0-th droplet only (go = 1): 

In this case the vector field d/dT^""^ is represented as 

d d 



_ + y]L(^(a)^ + 'ibJa)~] (31) 



(q2) Water is sucked from the a-th droplet and injected into the 0-th one with the unit 
rate: qp = —Sa/s, with no pumps in oil: 

To — Tk — 0, Sfs — —Sap 
d d 



9T(°) dS, 



(32) 



a 



In (29), (31) a/9T(") is to be understood not partial derivative but vector field 
in the space of g + 1 contours. By construction, it is an invariant vector field, i.e., it does not 
depend on the particular basis of harmonic functions and corresponding local coordinates. 
The same is true for d/dT^°'\ For a general process we have: 



Ql = Q^gf-^+T.Qj gf(^ +12P^-g^ (p-representation) 
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^ = <5oo^ + E + E (^-representation) 

The vector fields and relations between them are to be understood as acting on any physical 
quantity depending on the shape of the growing domain. It is important to stress that the 
linear superposition works only for processes of the same type (i.e., either of the p-type or 
g-typc). 

At fixed positions of oil sinks ILG spans a finite-dimensional subspace of an infinite 
dimensional variety of -|- 1- domains. In this subspace d/dT^°'^ and d/dT^"^ act as partial 
derivatives. For example, fix N points aj and consider the variety of contours which can be 
obtained from some initial configuration of droplets as a result of an ILG process with oil 
pumps at the points aj. The resulting shape of the droplets is uniquely determined (if no 
singularity occurs) by total amounts of oil sucked out from each point. This configuration 
space is A^-dimensional, and T^""^^ are local coordinates in it. Similarly, one may consider a 
more general configuration space, where additional parameters are amounts of water injected 
into each droplet. 



3 Analytic and algebro-geometric objects associated to 
the ILG 

In this section we describe analytic and algebro-geometric objects [24] which emerge in a 
description of an evolution of multiply-connected domains. 

3.1 Green function, harmonic measures and period matrix 

Pressure in the oil domain is expressed in terms of the following objects: 

• G{z,z'): The Green function of the Dirichlet boundary problem in Dqu. The function 
G{z,z') is symmetric and harmonic everywhere in Doii in both arguments except z — z' 
where G{z, z') = log \z — z'\ -\- . . . as z z'; besides, G{z, z') = if any of the variables 
z, z' belongs to the boundary. The Green function obeys the equation AG(^, z') — 
2nS^^\z-z'). 

• uja{z): The harmonic measure of the a-th boundary component. The function uja{z) is 
the harmonic function in Doii such that it is equal to 1 on dDa and vanishes on the 
other boundary curves. Thus the harmonic measure is the solution to the particular 
Dirichlet problem. The solution is given by 

<^a(^) = -77-/ dnG{z,C)ds, a = 0,l,...,g (33) 

SO the harmonic measure is the period of the differential d^Gdz. Obviously, the sum 
of the harmonic measures of all boundary components, dD^, which we call cycles, is 
equal to 1. In what follows we consider the linear independent functions u!a{z) with 
a=l,...,g. 
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• ^ap- The period matrix. Taking integrals of u!a{z) over nontrivial cycles, we define 



^aP^-^i dnUJaiOds, «, /3 = 1, . . . , (34) 
ZTT JdDp 

The matrix Q^iS is known to be symmetric, non-degenerate and positively-definite. It 
is called the period matrix because of its direct relation to periods of holomorphic 
differentials on the Schottky double of the domain Dou (see below). 

We also need the following "modified" objects, which are dual, with respect to the choice 
of the basis of canonical cycles (see Sec. 3.2) to the ones introduced above. 

• G{z,z'): The modified Green function [25] defined by 

G{z,z') = G{z,z') - ^ uja{z){n-%f,ujfs{z') (35) 

a,/?=l 

This function obeys the same equation AG{z,z') = 2n6^'^\z — z') and integrals of 
dnG over all the cycles 9Di, . . . , dDg are zero. However, instead of being zero on the 
boundaries, G takes there different constant values. 

• Lija{z): The modified harmonic measure is defined by 

ilj^iz) = -2j2{^-')apuJi,{z) (36) 

/3=1 

This is simply a linear combination of uJaS with domain- dependent coefficients such 
that ^ 

— i dnOJpds = 2Sap 

2tT JdDac 



The pressure field. Let us demonstrate how the solution for pressure in Doo is written 
in terms of the objects just introduced. For simplicity we do this assuming only one sink of 
oil with the power Q located at a point a {a — oo is also possible). 

The general solution for the pressure field p — p{z) with p — p^ on the boundaries reads 

9 9 

p{z) = 2QG{z, a) + Y, PaUJa{z) = 2QG{z, a) - QaOJaiz) (37) 

a=l a=l 

It is important to note that the rates qa completely determine pressures Pa in the water 
droplets and vice versa. Indeed, plugging (37) into (3), we have the relation 

2Q dnG {a, z)ds + ^Pp j^^ dnUJfj{z)ds = -iirqa 
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Using (33) and (34), we can write it either as a system of linear equations for (in the 
g-representation) , 

1 ^ 

X] = - Q'^a(a) (38) 



^^1 

or a system of hnear equations for (in the p-representation) , 

2j2i^)a}Qp-Pa-Qojaia) (39) 

13=1 

Since the g x g matrix Qai3 is non-degenerate, the system has a unique solution which is read 
from the equivalent "dual" system. Whichever the physical conditions in the water droplets 
are, pressure is given by 

9 9 

p{z) = 2QG{a,z)-2Q £ u;a{a){n-%pujp{z) + 2 £ qa{^-%pM^) (40) 

a,/3=l a,l3=l 

Specifying this formula for the elementary processes, we have: 

_ r 2G{a,z) for (pi) _ f 2G{a,z) for (ql) . . 

P^^^-\uj^{z) for(p2) ' P^^^-\u;a{z) for (q2) ^^'^ 

Variational formulas. Variations of the Green function and harmonic measures under 
infinitesimal deformations of the domain are described by remarkable formulas going back 
to Hadamard [26, 24]. Let Sn{^) be the normal displacement (with sign) of the boundary 
under the deformation counted along the normal vector at the boundary point ^, with the 
normal vector looking inside Doii, see Fig. 4. The variational formulas are: 

6G{z, z') = ^i d^z, dnG{z\ HO ds (42) 
27r JdD^ii 

SiUaiz) = ^ / dnG{z, dnUJc^iO Sn{0 ds (43) 
27r JdDM 

5^aP dnUJaiO dnMO HO ds (44) 

Small variations of the modified objects (with tilde) are described by exactly the same 
formulas (42), (43), where one should put tilde everywhere. 

These variational formulas arc easy to understand. Let us explain them on the example 
of the Green function. Since the Green function G{z,C,) vanishes if ^ belongs to the old 
boundary, the variation 6G{z, ^) is equal to the value of the new Green function on the old 
boundary, i.e., 5G{z,^) — —Sn{^)dnG{z,^), in the leading order. Now notice that SG{z,^) 
is a harmonic function (the logarithmic singularity cancels since it is the same for both old 
and new functions) with the boundary value —Sn{^)dnG{z,^). This function can be restored 
from its boundary value by solving the Dirichlet boundary problem. The r.h.s. of eq. (42) 
gives the result. 
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Figure 4: The normal displacement of the boundary. 

When the domain evolves with time, the Green function and harmonic measures become 
time dependent. The variational formulas allow one to find time derivatives of these functions 
given a local law of motion of the boundary. This is the way how we derive partial differential 
equations for the ILG below. 

3.2 The Schottky double 

The growing domain Doii is a Riemann surface with a border. From mathematical point of 
view, it is more convenient to work with compact Riemann surfaces without border than 
with bordered domains. Given a planar domain with holes, like Don, endowed with the 
holomorphic coordinate z, it may be thought of as a "half" of a closed Riemann surface. 
Another half, an antiholomorphic "copy" of Doii with coordinate z, is glued to the first copy 
along the boundaries dDa- Besides, each copy of Doii should be compactified by adding a 
point at infinity. The resulting compact Riemann surface without boundary is called the 
Schottky double, or simply the double of the planar bordered domain (sec, e.g. [24]). 

More precisely, the Schottly double of a bordered surface, is a compact Riemann surface 
without boundary endowed with an antiholomorpic involution such that the boundary of 
the initial domain is the set of fixed points of the involution. The Schottky double of the 
□oil can be thought of as two copies of it ("upper" and "lower" sheets of the double) glued 
along the boundaries Uf^^^dDa, with points at infinity added (cxo and ob). The holomorphic 
coordinate on the upper sheet is z, the same as in Don, while the holomorphic coordinate on 
the lower sheet is z. The genus of the double is equal to the number of water droplets minus 
1. The notion of the double was implicitly used by Richardson in his studies of viscous flows 
in multiply-connected fluid regions [7, 8]. 

A meromorphic function on the double is a pair of meromorphic functions /, / in Doq 
such that f{z) = f(z) on the boundary. Similarly, a meromorphic differential on the double 
is a pair of meromorphic differentials f{z)dz and f{z)dz such that f{z)dz = f{z)dz along 
the boundary curves. The Schwarz reflection principle says that any meromorphic differen- 
tial di'{z) on the upper sheet such that it is purely imaginary along the boundary, can be 
meromorphically extended to the lower sheet as —dv{z), so for each pole of such a globally 
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Figure 5: The two halfs of the Schottky double (g — 2). 

defined differential there is a "mirror" pole on the opposite sheet. 

To proceed, one has to choose a basis of a- and b-cycles on the double having the canonical 
intersection form o = o = 0, Sa obp = Sa/3- In general, for an abstract Riemann 
surface, there is no preferred choice of the basis. However, when the surface is the double of 
a planar domain, like in our case, we may fix the following two distinguished ( "dual" ) bases. 

• The b-cycles are just boundaries of the holes bo, = —dD^ for a — 1,. . .,g. (Note, 
however, the negative clockwise orientation.) The a^-cycle connects the boundary of 
the «-th water droplet with the 0-th one. To be more precise, fix points on the 
boundaries, then the a„-cycle starts from ^o- goes to on the "upper" sheet and 
returns following the same way on the "lower" sheet. 

• In the "dual" basis, we interchange a- and b-cycles: aa = — ba, b^ = a^. The minus 
sign is necessary to preserve the anti-symmetric intersection form. 

These two choices of the basic cycles correspond to the LG dynamics in the p and q- 
represent ations . 

3.3 Differentials on the Schottky double 

Holomorphic differentials. Having fixed the basis of cycles, we can introduce normalized 
holomorphic differentials (abelian differentials of the first kind) as differentials of holomorphic 
parts of the harmonic measures. Let us represent uJa as the real part of a holomorphic 
function: = Wa{z) + Wa{z), where Wa{z) are holomorphic (in general multivalued) 
functions in Doii. The differentials dWa — dzCUadz are then holomorphic in Doii and purely 
imaginary along all boundary contours. So they can be extended holomorphically to the lower 
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p -representation 



q - representation 



Figure 6: Canonical cycles on the Schottky double in p and q-representations. 



sheet as —dWa{z). They form the canonically normahzed basis of holomorphic differentials 
on the double w.r.t. the a-cycles: 

/ dWp = j^" dWp{z) + i-dW^) = 2Re j^" dWp{z) = doop = 

The matrix of b-periods of these differentials reads (cf. (34)) 

Tap = f dWp = f dnUJ/sds = iirQaiB (45) 

Similarly, the differentials dWa = d^u^adz, constructed in the same way from the modified 
harmonic measures, are normalized w.r.t. the a-cycles^: 

j> dWfj = 2TCi 5af3 

The b-periods of the differentials dWa are: 

I dWp = -2{n-')af3 

J ha 

Meromorphic differentials [24]. The simplest meromorphic differential on a Riemann 
surface is the differential having only two simple poles with residues ±1 (the abelian differ- 
ential of the third kind). We will especially need the differentials whose only simple poles 
are at the mirror points of the double. They are naturally associated with the oil pumps. 

^Normalization of the holomorphic differentials differs by 27rz . One might work using the unified normal- 
ization in both p and g-representations but in that case, as a price for the unified notation, some artificial 
imaginary units enter formulas for physical quantities. 
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Specifically, consider the differential dW'^"''"'\z) — 2dzG{a, z)dz on the upper sheet. Prom 
the properties of the Green function it follows that this differential has the only simple pole 
with residue +1 at the point a. Since along all the boundaries 

= d,G{a, z) = d,G{a, z)dz + d^G{a, z)dz = 2Re(diy("'")(z)) = 

this differential can be meromorphically extended to the lower sheet as —2dzG{a, z)dz, and 
so it has a simple pole with residue —1 at the mirror point d on the lower sheet. In particular, 
the differential dW'^°°'°°\z) — 2dzG{oo,z)dz has simple poles at the two infinite points of 
the Schottky double. The differentials dW'^"''"'^ arc constructed in the same way out of 
the modified Green function. Note that the so defined differential dH^^"'"-* are canonically 
normalized: 



means Re / dzG{a, z)dz — (p-representation) 
j> dW^'-"'"^ = that means j> dzG{a, z)dz — (g-representation) 

Ja^ JdD„ 



The abelian differentials of the second kind, dW^^\ have only one pole of order A; + 1 at 
the point b. These differentials can be explicitly defined by expansion of the Green function 
or the differential dW^""'^^ in a Taylor series in a — 6 in the vicinity of the point b. This 
differential is normalized, i.e., dW^'^^ — but to define it uniquely one has to fix the 
principal part at the pole. 

Let us introduce the differentials dW^°°^ with the pole at infinity. Expanding the Green 
function around oo on the upper sheet, and using the basis of analytic functions ipk (15), we 
write 

2dzG{a, z)dz = 29,^(00, z)dz - ^ ^ [ijk{a)dwt\z) + j;jajdwt\z)) (46) 

fc>i 

Here dW^^^ is the normalized differential holomorphic everywhere except infinity (on the 
upper sheet) with the principal part 

dwi°°\z) = kdxk{z) + Oiz-"^) as z ^ oo 

On the lower sheet, this differential is regular. Similarly, dwl^'^ has a pole at infinity on the 
lower sheet of the double and is regular on the upper sheet. The differentials are defined by 
the following contour integrals: 

kdz f 

dz f Xk{Odc^G{z, ()d( (on the upper sheet) 

(47) 

kdz f 

— r dz f Xk{C)9(;G{z, ()d( (on the lower sheet) 

The contour here encircles oo but not the point z. The differential dwl°°^ is defined as the 
complex conjugated differential {—dwjf^'') on the opposite sheet of the double. 

The differentials introduced above are summarized in the table: 



dwt^ 
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Type 


Notation 


Explicit form on the up- 
per/lower sheet 


1 

Jaf3 




First kind 




dzUJa{z)dz 
— dzUJa{z)dz 






Second kind 




see (47) 
see (47) 





kf dzXk{z)dzUJa{z) 


Third kind 




2dzG{a, z)dz 
— 2dzG{a, z)dz 








Correspondence between times and differentials. Each elementary flow is naturally 
"coupled" with a meromorphic or holomorphic differential on the Schottky double. This 

differential can be defined in terms of pressure as dzp{z)dz. Equivalently it is a unique 
analytic continuation of the differential dtz{a,t)dz — dtz{a,t)dz to Dqu from the boundary 
(here a is any parameter on the contours). The equivalence of the two definitions follows 
from the Laplacian growth equation (14) (Sec. 2.2). Explicitly, the coupling is: 

I , , r(") I > dWa , Tk I > dWl'^'^ 

In particular, 7"^°°) = Tq + const, corresponding to an oil sink at infinity, is coupled with 
the differential dW^°°'°°\ The abehan integral 

W{z)^W^'^'^\z)= r dW^°^^'^^ (48) 

-'Co 

where Co is an arbitrary fixed point in Dqh, has a constant real part along each component 
of the boundary and W{z) — log 2; + 0(1) as 2; — > 00. For the process (ql) with the sink at 
infinity the function w{z) = e^^^^ conformally maps the domain Dqu onto the exterior of a 
disk of radius tq = \w{C,o)\ {^0 is any point on bo) with g concentric arcs removed. Radii and 
endpoints of the arcs depend on time. The radius, Va, is related to pressure in the a-water 
droplet by Pa = 21og(ra/ro). 

4 Equivalence of the ILG and Whit ham dynamics 

We describe an evolution of the Riemann surface through a time dependence of the canonical 
abelian differentials and abelian integrals on the Schottky double. One may realize them 
as holomorphic or meromorphic functions of the coordinate z on the physical plane. These 
functions encode the shape of the interface at any instant of time. In this respect they serve 
as substitutes for the time-dependent conformal map to the unit disk in case of a single 
droplet. 
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We represent the dynamical equations of the ILG as relations between variations of the 
abelian integrals under simultaneous action of different oil and water pumps. Treating T^"^ 
and T(") as an independent "time" variables the evolution takes the form of a hierarchy of 
partial differential equations in many variables. We recognize them as a universal Whitham 
hierarchy [14, 15] in the (extended) moduli space of genus g Riemann surfaces (see Sec. 4.3). 



4.1 Exchange relations 

Let us calculate the time derivative of the Green function G{a, b) in the elementary process 
with the oil pump at some third point c. Using the Hadamard formula (42), and the fact 
that for this process Sn{^) — —^dnG{c,^)5T^'^\ we have: 



dG{a, b) __1 r 



dr,G{a,OdnG{b,OdnG{c,Ods 



Remarkably, the result is symmetric under all permutations of the points a, b, c. For a single 
connected domain this equation has been obtained in Ref. [27]. In a similar way, we find, 
using (43), (44): 



dnUJa{09nG{b, ^)dnG{c, ^)ds 



and 



47r /aD„ii 

as well as analogous formulas for derivatives w.r.t. T^"\ Observing the symmetry of the 
right hand sides of these equations, one may write them as local relations connecting the 
time derivatives of the Green function and harmonic measures. 

In this way we obtain the following fundamental exchange relations: 

dG{a, b) dG{b, c) dG{c, a) 



dTi") 5T(") dT(^^ 
duja{a) dojaib) 



(49) 

(50) 
(51) 



which encode dynamics of the Laplacian growth with zero surface tension. There is also the 
relation which connects derivatives of the Green function and harmonic measure: 

n ^G{a,b) _ dujaja) 

In the same way, one may extend this list to include derivatives of the period matrix: 

^dup{a) dVLaB dQafs dVLg^ dQ^a 



grpia) Qrjn(a) ' Qrjn^-f) Qrp{a) Qrp{l3) 

For the dual objects, the relations (49) - (52) remain the same. 



(53) 
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4.2 The Whitham equations 



Whitham equations are holomorphic parts of the exchange relations. Taking, for instance, 
dz of both sides of the exchange relation dj^(b)G{a, z) = drp(a)G{b, z), we get dj'(b)dW^°''"'\z) — 
dj'(a)dW'^^'^\z) which is an equation of the Whitham hierarchy. The full list of Whitham 
equations obtained in this way reads: 

^ -dW^'f\z) ^ J-^,dW^^'-^\z) (54) 



^ -dW^'^'-^^z) = ^dW^{z) (55) 



^dWpiz) = ^dWaiz) (56) 



The derivatives are taken at constant z. The hst can be further enlarged by adding the 
equations containing T^-derivatives: 

^ -dWt\z) = 4^dWi-\z) 



and so on. They can be obtained from the generating equations (54) - (56) by expending 
them in a series around infinity. 

We note that there are a few other equivalent ways to write the Whitham equations 
[15]. Let us present a more invariant formulation, which does not rely on the choice of the 
distinguished local coordinate z in Doii. The Whitham equations in the invariant form are 
naturally written in an extended "moduli space" of contours. Points of this extended space 
are sets of data of the form 

(g + 1 nonintersecting boundary contours; a point in Doii^ 

Local coordinates in this space are parameters T$ or TS and a local coordinate A in the 
domain Dou. 

The local coordinate A may depend on the shape of the domain. In particular, one may 

choose A to be one of the abehan integrals, say W{z) = W^°°'°°\z) with dW — (48). 
Then the function w{z) = e^^^^ is single- valued in Doii- It is a good local coordinate in Doii 
everywhere except for the points where dW^°°'°°^ = 0. It can be shown that all these points 
belong to the boundaries bi, . . . , bg and there are exactly two such points on each boundary. 
Under the conformal map w{z) these points are taken to the endpoints of the concentric arcs 
(see the end of Sec. 3.3). Treating all other abelian integrals as functions of W rather than 
z^z{W,T), 

J Cn 



'Co 

we represent the Whitham equations (54) - (56) in the form 
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where A stands for a or a and 



dW dTo dW dTo 



has the form of "Poisson brackets". These equations are consistency conditions for the 
system of evolution equations 

^ = {^'". (58) 
Eqs. (57,58) constitute the Whitham universal hierarchy of the sohton theory. 



4.3 Whitham equations in sohton theory 

Whitham equations as modulation equations. Integrable partial differential and dif- 
ference equations of soliton theory are known to possess a rich family of periodic exact 
solutions depending on continuous parameters. To be more definite, we start our discussion 
with a (1 + 1) -dimensional integrable evolution equation of the form Ut = P{u, u^, ■ ■ ■) (e.g., 
the KdV equation) . Exact periodic solutions have the form 

u{x, t) = Uo{\Jx + Yt + Z\I) 

where U, V, Z are g'-dimensional constant vectors with components {C/q} = (C/i, . . . , Ug), 
etc, and uq{Z) is a periodic function of any component Z^. This function, and all the 
vectors, depend on the set of parameters / = (/i, . . . ,/m)- Each periodic solution can be 
constructed starting from a Riemann surface. / stands for the moduli of a Riemann surface 
and the vectors U, V are b-periods of certain normalized meromorphic differentials, dW^^'^ 
and dW^^\ on the Riemann surface, with prescribed singularities at infinity [1]. 

In a number of physical problems one is interested in slowly modulated periodic solution, 
rather than just periodic. A nonlinear WKB method or the Whitham averaging method 
allows one to construct more general solutions of the same integrable equation using the 
function uq as a leading term of the asymptotic expansion 

u{x, t) = uo[e-^S{X, T) + Z(X, T) \ I{X, T)) + eui{x, t) + e'^U2ix, t) + . . . 

where £ is a small parameter and the parameters / now depend on the slow variables X — ex, 
T — et. The original variables x, t are called fast variables. If the vector-valued function S 
obeys the equations 

dxS = V{I{X, T)) , drS = V{I{X, T)) 

then the leading term agrees with the original solution up to first order in e. All the higher 

corrections can be found by solving non-homogeneous linear equations whose homogeneous 
part is the original equation linearized on the background of the exact solution uq [12]. 

We see that the so constructed solution u{x,t) describes the original fast periodic oscil- 
lations, modulated, on a larger scale, by a slow drift in the space of exact periodic solutions. 
The equation, which describes the drift I{X,T) are called Whitham equations. For the 
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particular example discussed above the equation to determine I{X,T), written in a proper 
local parameter reads dxdW^*^ — dxdW^^h This form of the Whitham equations was first 
observed in [13] for the KdV hierarchy. Implicitly, through the dependence of the canonically 
normalized differentials on the slow variables, they describe the drift in the moduli space of 
Riemann surfaces and thus the dependence I{X, T). The Hamiltonian approach to Whitham 
equations for (1 + 1) -dimensional systems was developed in [28]. A universal Whitham hi- 
erarchy in a general setting of multi-dimensional integrable equations was suggested in ref. 
[14] . An invariant formulation of the Whitham hierarchy, independent on the choice of local 
coordinates, was given in [16]. 

The idea of [14] was to obtain equations describing the slow drift in the space of exact 
solutions from the condition that next-to-leading terms of the asymptotic series u{x,t) be 
uniformly bounded on large scales. In general, the asymptotic series becomes unrehable on 
scales of order e~^, i.e., the corrections become large. The main result of [14] is that the 
Whitham equations follow if one requires that just the next term of the series, ui{x,t), be 
uniformly bounded for all x,t. 



The universal Whitham hierarchy. In a more general multi-dimensional hierarchy of 
soliton equations (like the KP hierarchy or the 2D Toda lattice hierarchy, or their difference 
counterparts), one has a family of "times" Ia and a family of "potentials" u^^\ evolving 
with times tA- (Here A and N belong to a case-dependent, generally infinite set of indices.) 

Exact periodic solutions of the hierarchy are constructed from a given time-independent 
Riemann surface F with some additional data on it. With each time tA one associates a 
meromorphic differential, dWA, on F normalized with respect to, say, a-cycles: fg^ dWA — 0. 
Let U^^) be the vector of b-periods of this differential: 

C/(^) = / dWA 

Then the exact solution has the form {u is one of it*^-^)'s): 

u{{tA}) = uo(EU^^V)^A + Z(^) 1 1) + co{I) 

A 

where Uq is a certain oscillating periodic function (the second logarithmic derivative of the 
Riemann ^-function) and cq is a constant. Like in the previous example, one may try to con- 
struct a more general oscillating, but not periodic solutions with slowly varying parameters 

u = uo{e-'S{{TA}) + Z{{Ta}) \ I{{Ta})) + co(/({T^})) + eui{{tA}) + ... (59) 

Here Ta — sIa are slow times and S is a vector function such that 9t^S = U(^)(/({T4})). 
The uniform boundness of the first correction ui for all times imphes the hierarchy of 
Whitham equations [15], 

OWa ( OWb _ dWc\ dWB ( dWc _ dWA dWc ( dWA _ dWA ^ , . 

dX [ dTc OTb dX [ OTa dTc dX [ OTb OTa J ^ ' 
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valid for all possible values of the indices A,B,C. Here A is any local parameter, all the 
abelian integrals being regarded as functions of A. Choosing one of the indices, say C, to be 
(oo, 6b) and A = 1^ = w'^°°'°°\ one gets the Whitham equations in the form (57). 

While averaging the solution (59) over fast oscillations, {uq) vanishes 

(«)(ra) = co(ra) 

In the context of Laplacian growth, T4 is T*^"-* or T^°'\ and cq is the Green function G{a, b). 
Here the points (a, b) label the potential. Thus the Laplacian growth can be thought of as 
a physical realization of the slow drift in the moduli space of Riemann surfaces. 

5 Special classes of solutions of the Whitham hierarchy 
and Laplacian growth of algebraic domains 

In this section we briefly discuss particularly important families of growing domains. They 
correspond to special solutions to the Whitham hierarchy called "algebraic orbits" [15]. 

5.1 Algebraic and abelian domains 

Let us recall that the time derivative of the Cauchy transform of the oil domain (10) is a 
globally defined rational function. Consider the class of domains whose Cauchy transform is 
a single globally defined meromorphic (i.e., rational) function in the plane. In other words, 
each function defined originally in the domain is extendable to a single rational 
function defined everywhere in the plane, same for different a. In this case one may forget 
about the index a and deal with the single function h{z). Eq. (10) tells us that if the initial 
fluid region is from this class, then it remains to be in this class in the process of the LG 
evolution. The evolution may only add new poles or change residues of the existing ones. As 
is pointed out in [8], in order to prepare such an initial condition, one may inject oil through 
points into a cell initially filled by water. 

The domains whose Cauchy transform is a globally defined rational function are called al- 
gebraic [11] or quadrature domains [29]. Some illustrative examples in the multiply-connected 
case can be found in [8] and [30] . This class appears to be quite representative and important 
since any domain with smooth boundary components can be approximated by quadrature 
domains (see [31] for the proof). In the simply-connected case, the quadrature domains are 
images of the unit disk under conformal maps given by rational functions. Their time evo- 
lution is described by rational solutions of the Laplacian growth equation. These solutions 
are sometimes known to develop cusp-hke singularities within finite time [19]. 

A more general class of domains can be defined by imposing the above condition not 
on the /i+ itself but on its ^-derivative. Namely, suppose that each differential dh'^{z) is 
extendable to a meromorphic differential in the plane, and they coincide for different a's. In 
this case itself may be a multivalued analytic function with logarithmic branch points. In 
[11], such domains were called abelian domains. They can be produced from the quadrature 
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domains by the oil sucking from linearly extended sinks. In the simply- connected case, their 
evolution is described by logarithmic solutions of the Laplacian growth equation [32, 33]. 

5.2 The Schwarz function 

For domains with analytic boundaries, and for algebraic domains in particular, the Cauchy 
transform allows one to introduce the Schwarz function of the boundary contours, which 
proved to be very useful for analyzing the LG dynamics in the simply-connected case [2]. 
Given a closed contour on the plane, the Schwarz function [29, 34] is defined as the analytic 
continuation of the function z away from the contour. Let us denote it by S{z). According 
to the definition, S{z) is a function analytic in some neighborhood of the curve such that 

S{z) — z on the curve (61) 

From the continuity of the Cauchy transform we have z = h^{z) — h^{z) on the a-th 
boundary, so S^"\z) = h^{z) — h~{z) is the Schwarz function of the a-th boundary curve. 
In general all these functions are different. However, as it directly follows from the definition, 
all the boundary contours of algebraic or abelian domains have a common Schwarz function, 
S{z) = S^°'\z) for any a, and the differential dS{z) is meromorphic in Don, i.e., it has there 
only a finite number of isolated poles. 

In the case of algebraic or abelian domains, one may decompose the Schwarz function into 
the sum S{z) = S^{z) +S-{z), where the function = h'^{z) (for any a) is analytic inside 
water droplets while S'_ = —h~{z) is analytic in D^n and vanishes at infinity. Combining the 
time derivatives of and h~ (see (8) and (9)), we find the time derivative S{z) — dtS{z;t) 
of the Schwarz function at constant z: 

S{z) = h+{z) - h-{z) = -d,p{z) (for z in Don) (62) 

where by ha{z) we mean the analytic continuation of this function to Dqu. 

In D^s, the analytically continued function S{z) has more complicated singularities. In 
case of general position, they are branch points r]i of order two (with cuts between them). 
One may think of the algebraic curve underlying the solution to the LG problem as the 
Riemann surface of the function S{z). The Whitham equations can then be equivalently 
represented in the form of equations for the endpoints of the cuts, which describe their time 
dependence rji — r]i{T). The equations are: 

% ^ ( dW(^\z) \ dru 
dTi^) [ dWiz) J^^^^dTo ^^""^ 

where the coefficients in front of drji/dTo in the r.h.s. are expressible as (in general qiiite 
complicated) functions of the ?7fc's. Basically these equations mean that the Schwarz function 
takes finite values at the branch points r]i. 
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Figure 7: The integration path from = to z^- 



5.3 The generating differential 

For the class of algebraic domains, one may introduce a distinguished meromorphic differ- 
ential on the Schottky double. Recall that the Schwarz function S{z) is a meromorphic 
function in Doii, if Doii is an algebraic domain. Therefore, one may treat S{z) as a function 
on the Schottky double extending it to the lower sheet as z. In this case the differential 
dS — S{z)dz is extendable to a meromorphic differential on the double. Its explicit form on 
the lower sheet is zdS{z). 

The differential dS plays an important role in the theory. We call it the generating 
differential because it generates the complete set of local coordinates in the space of algebraic 
multiply-connected domains. Indeed, the harmonic moments are 



The areas of water droplets are periods of this differential over b-cycles: 

Sa^-^J, dS (65) 



2m Jh, 

The "dual" variables, are expressed through periods of the same differential over the 
a-cycles (or, what is the same, over b-cycles). To show this, we write: 

^a= d(p{z) ^2Re h{z)dz^ 
Jo Jo 

= 2Re h^{z)dz + J'" h'^{z)dz + J^" {z + h' {z))dzj 



where the integral from to Za goes along a path containing the two points and ^o, ia are 
intersection points of the path with the boundary curves bo, ba (see Fig. 7). Adding and 
subtracting J^^ h'^[z)dz (which is well-defined for algebraic domains), we rewrite this in the 
form ^ 

= 2Re / " S+{z)dz - n„ 

^0 

where 

Ha = 2Re / \S{z) - z)dz = / " (S{z)dz - zdS{^) ^ i dS (66) 

Jin Jin ^ ' J^c, 
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Note that does not depend on the points ^o, Ca- Finally, the times T^"^ associated with 
a sink of oil at the point a are expressed as 

= / ds (67) 

2m Jca 

(Here is a small contour encircling the point a.) 

The table of differentials from Sec. 3.3 can be continued by including the generating 
differential: 



lype 


Notation 


Explicit form on the up- 
per/lower sheet 


£ 




Generating 
differential 


dS 


S{z)dz 




—2'KiSa 






zdS{z) 







Moreover, partial derivatives of the generating differential w.r.t. the times T^°'\ T^"'^ 
coincide with the canonical meromorphic differentials: 

dz = -dW^'^'^^z) , = -dW^(z) (68) 

This follows from (62) after substituting the pressure filed for the elementary processes in 
terms of the Green function and harmonic measures (see (41)). For algebraic domains, the 
Whitham equations follow from the existence of the generating differential and (68). 

The generating differential can be represented as 

dS^dk-Y^T^'^MW^^^ (69) 

A 

where rfA is a differential with T'^^-'-independent singularities. For algebraic orbits, it is 
a fixed meromorphic differential (possibly with time-independent jumps). In more general 
cases dK has more complicated analytic properties. Presumably, it can be defined as a 
solution to a 9-problem. The expansion (69) and equation dS{\^ = for all zeros As of the 
differential dz on the lower sheet of the Schottky double, where it has the form dS{z), arc key 
relations which imply (68). Indeed, from the latter condition it follows that the differential 
dAdS has no singularities at the points A^. Then (69) implies that this differential has the 
same singularities and periods as —dW^'^^ Hence, they do coincide. 

General algebraic orbits of Whitham equations for higher genus Riemann surfaces in the 
sense of [15] correspond to the case when the z-derivative of the Schwarz function extends 
to a meromorphic function on the double (equivalently, when the differential dS{z) extends 
to a meromorphic differential on the double). 
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6 Conclusion 



In short, the main message of this work is that the variables in which the Laplacian growth 
with zero surface tension becomes hnear, for arbitrary connectivity of the growing domain 
and arbitrary configuration of pumps, are the Whitham "times" defined in [15]. The latter 
are special local coordinates on the extended moduli space of Riemann surfaces. Conservation 
or linear dependence on time of harmonic moments of the growing domain, known before as 
a characteristic feature of the idealized Laplacian growth, is a particular case of this result. 

The Whitham equations are partial differential equations for canonical holomorphic and 
meromorphic differentials on Riemann surfaces regarded as functions of the local coordinates 
in the moduli space. Solutions to the Whitham equations allow one to find the differentials 
and abelian integrals as functions of time and reconstruct dynamics of the interface. 

The Whitham equations are often regarded as integrable ones, though not in the Liou- 
ville sense. When speaking about integr ability of Whitham equations, one means mainly a 
possibility to actually integrate them by representing a solution in the form of an implicit 
function of independent variables (the hodograph method). 

The Whitham equations appear in soliton theory in different contexts. First, when one 
looks for solutions of soliton equations other than periodic. At some regimes, these solutions 
are well approximated, on small space-time scales, by the periodic exact solutions of the 
algebro-geometric type. When fast oscillations of the periodic solutions are averaged or 
smoothed out, Whitham equations appear as modulation equations written for moduli of 
the Riemann surface parametrizing the algebro-geometric periodic solutions. An important 
special case of Whitham equations appears if one neglects dispersion in nonlinear soliton 
equations. The latter case provides the most direct link to the Laplacian growth of simply- 
connected domains. This link was explored in Ref. [2]. What is perhaps the most important 
conclusion, the Whitham equations describe a proper evolution of the Riemann surface built 
upon a growing interface. 

The relation between the growth problem and modulated periodic solutions to soliton 
equations is two-fold. The Laplacian growth may serve as a simple illustrative physical model 
of the Whitham dynamics of complex curves. Vice versa, the methods developed in soliton 
theory may help to understand growth in a singular (turbulent) regime, i.e., in a vicinity of 
cusp formation or coalescence and break-up of droplets, providing an effective account of the 
surface tension effects near singular points. 
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